# BASH Interface packages
import os
import glob
import subprocess
# Dataframe pacakegs
import pandas as pd
import numpy as np
# Jupyter plugins
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import warnings
import math
# Plotting packages
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
warnings.filterwarnings("ignore")
# Theme stuff
# from jupyterthemes import jtplot
# jtplot.style(theme='grade3')
# jtplot.style(context='talk', fscale=1.4, spines=False, gridlines='--')
# jtplot.style(ticks=True, grid=False, figsize=(6, 4.5))
# Defines a few functions that:
# Filters through the pd dataset
def filter_class(basename, ionname, energy):
return(list(data[data["base_name"] == basename][data["ion_class"] == ionname][energy]))
def filter_class2(basename, ionname, energy):
return(list(homodata[homodata["base_name"] == basename][homodata["ion_class"] == ionname][energy]))
# Runs bash commands
def runbash(command):
proc = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, check=False, text=True, executable='/bin/bash')
output = proc.stdout.strip('\n')
if output != "":
return(output)
def Eh2eV(df, column):
for icolumn in column:
for i in range(len(df[icolumn])):
df[icolumn][i] = df[icolumn][i] * 27.2116
def filter_class_pairs(basename, ionname, energy):
return(list(datapairs[datapairs["base_name"] == basename][datapairs["ion_name"] == ionname][energy]))
def filter_class2_pairs(basename, ionname, energy):
return(list(homodatapairs[homodatapairs["base_name"] == basename][homodatapairs["ion_name"] == ionname][energy]))
def filter_ions_pairs(ionname, energy):
return(list(datapairs[datapairs["ion_name"] == ionname][energy]))
# The lists to cycle through
basenames = ["ethane", "ethene", "ethyne", "benzene", "furan", "pyrrole"]
ionnames = ["BF2(CF3)2", "BF4", "C(CN)3", "N(SOOCF3)2", "Cl" ,"N(NC)2", "C1MPyr", "ImidMe2", "NMe4", "PMe4", "Na", "Li", "Pd", "Ni", "tempo"]
classnames = ["cation", "anion", "metal"]
classnamespairs = ["high", "low"]
energynames = ["electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "total", 'dhomo', 'dlumo', 'homo', 'lumo', 'gap', 'nbo_homo', 'nbo_lumo', 'nbo_gap', 'dplusel']
Since removing from the server it was running on with complptete access to live calcualtions as they're coming in, the data has been exported to csv files to be imported cleanly.
I could have converted these blocks to raw, but then my code hiding extension wouldn't hide them.
# homodata = pd.read_csv("./IonsData/homo-lumo-base.csv")
# # Sets the dataframe column type to float
# homodata.homo = homodata.homo.astype(float)
# homodata.lumo = homodata.lumo.astype(float)
# lone = pd.DataFrame(columns=['base_name', 'ion_class', 'ion_name', 'homo', 'lumo'])
# lone = lone.append(homodata[homodata["base_name"] == "alkyl"], ignore_index = True)
# lone = lone.append(homodata[homodata["base_name"] == "anions"], ignore_index = True)
# lone = lone.append(homodata[homodata["base_name"] == "cations"], ignore_index = True)
# lone = lone.append(homodata[homodata["base_name"] == "aromatic"], ignore_index = True)
# droplist = []
# for i in reversed(range(len(list(homodata.base_name)))):
# if list(homodata.base_name)[i] == "alkyl" or list(homodata.base_name)[i] == "aromatic" or list(homodata.base_name)[i] == "cations" or list(homodata.base_name)[i] == "anions":
# droplist = droplist + [i]
# homodata = homodata.drop(droplist)
# # Initalise the dataframe
# data = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'electrostatics', 'exchange', 'induction', 'dispersion', 'total', 'charge_transfer', 'dhomo', 'dlumo', 'bhomo', 'blumo', 'plotsize'])
# # This runs through the files and pulls out important data
# infiles = glob.glob("../4-SAPT/*/*.out")
# for i in infiles:
# if "slurm" not in i:
# fullfilepath = runbash("readlink -f '" + i + "'")
# filename = runbash("/bin/basename -s .out '" + fullfilepath + "'")
# ion_name = runbash("echo '" + filename + "' | cut -d'-' -f2")
# base_name = runbash("echo '" + filename + "' | cut -d'-' -f1")
# sort = {"ethane" :1, "ethene" : 2, "ethyne" : 3, "benzene" : 4, "furan" : 5, "pyrrole" : 6}
# classes = [["BF2(CF3)2", "BF4", "C(CN)3", "N(SOOCF3)2", "Cl" ,"N(NC)2", "tempo"],["C1MPyr", "ImidMe2", "NMe4", "PMe4", "Na", "Li"],["Pd", "Ni"]]
# if ion_name in classes[0]:
# ion_class = "anion"
# elif ion_name in classes[1]:
# ion_class = "cation"
# elif ion_name in classes[2]:
# ion_class = "metal"
# electrostatics = runbash("cat '" + i + "' | grep 'Electrostatics ' -m1 | tr -s ' ' | cut -d' ' -f7")
# exchange = runbash("cat '" + i + "' | grep 'Exchange ' -m1 | tr -s ' ' | cut -d' ' -f7")
# induction = runbash("cat '" + i + "' | grep 'Induction ' -m1 | tr -s ' ' | cut -d' ' -f7")
# dispersion = runbash("cat '" + i + "' | grep 'Dispersion ' -m1 | tr -s ' ' | cut -d' ' -f7")
# total = runbash("cat '" + i + "' | grep ' Total SAPT2+ ' -m1 | tr -s ' ' | cut -d' ' -f8")
# charge_transfer = runbash("cat '" + i + "' | grep 'SAPT Charge Transfer' | tail -n 1 | tr -s ' ' | cut -d' ' -f9")
# newdata = {"sort" : sort[base_name], "base_name" : base_name, "ion_class" : ion_class, "ion_name" : ion_name, "electrostatics" : electrostatics, "exchange" : exchange, "induction" : induction, "dispersion" : dispersion, "total" : total, "charge_transfer" : charge_transfer, "dhomo" : 0, "dlumo" : 0, "plotsize" : 0}
# data = data.append(newdata, ignore_index=True)
# # Drops any frames that have a "None" value in charge_trasfer (the last SAPT calculation)
# data = data[data['charge_transfer'].isnull() == False]
# # Sets the dataframe column type to float
# data.electrostatics = data.electrostatics.astype(float)
# data.exchange = data.exchange.astype(float)
# data.induction = data.induction.astype(float)
# data.dispersion = data.dispersion.astype(float)
# data.total = data.total.astype(float)
# data.charge_transfer = data.charge_transfer.astype(float)
# data.dhomo = data.dhomo.astype(float)
# data.dlumo = data.dlumo.astype(float)
# data.bhomo = data.bhomo.astype(float)
# data.blumo = data.blumo.astype(float)
# data.plotsize = data.plotsize.astype(int)
# data = pd.merge(data, homodata, on=['base_name', 'ion_name', 'ion_class'])
# warnings.filterwarnings("ignore")
# for i in range(len(data.base_name)):
# for j in range(len(lone.homo)):
# if data.base_name[i] == lone.ion_name[j]:
# data.dhomo[i] = data.homo[i] - lone.homo[j]
# data.dlumo[i] = data.lumo[i] - lone.lumo[j]
# data.bhomo[i] = lone.homo[j]
# data.blumo[i] = lone.lumo[j]
# for i in range(len(data.base_name)):
# if math.isnan(data.total[i]) == False:
# data.plotsize[i] = abs(int(data.total[i]))
# Eh2eV(data, ["lumo", "homo", "dhomo", "dlumo", "bhomo", "blumo"])
# Eh2eV(lone, ["lumo", "homo"])
# lonedata = pd.read_csv("./IonsData/homo-lumo-base.csv")
# lonepairs = pd.DataFrame(columns=['base_name', 'ion_class', 'ion_name', 'homo', 'lumo', 'gap', 'nbo_homo', 'nbo_lumo', 'nbo_gap'])
# lonepairs = lonepairs.append(lonedata[lonedata["base_name"] == "alkyl"], ignore_index = True)
# lonepairs = lonepairs.append(lonedata[lonedata["base_name"] == "aromatic"], ignore_index = True)
# # Sets the dataframe column type to float
# lonepairs.homo = lonepairs.homo.astype(float)
# lonepairs.lumo = lonepairs.lumo.astype(float)
# lonepairs.drop('ion_class', axis=1, inplace=True)
# for i in range(len(lonepairs)):
# lonepairs.base_name[i] = lonepairs.ion_name[i]
# lonepairs.ion_name[i] = "lone"
# homodatapairs = pd.read_csv("./IonsData/homo-lumo-pairs.csv")
# # Sets the dataframe column type to float
# homodatapairs.homo = homodatapairs.homo.astype(float)
# homodatapairs.lumo = homodatapairs.lumo.astype(float)
# homodatapairs.drop('ion_class', axis=1, inplace=True)
# # Initalise the dataframe
# datapairs = pd.DataFrame(columns=['sort', 'base_name', 'ion_name', 'electrostatics', 'exchange', 'induction', 'dispersion', 'total', 'charge_transfer', 'dhomo', 'dlumo', 'bhomo', 'blumo', 'plotsize'])
# # This runs through the files and pulls out important data
# infiles = glob.glob("../7-IonPairsSAPT/*.out")
# for i in infiles:
# if "slurm" not in i:
# fullfilepath = runbash("readlink -f '" + i + "'")
# filename = runbash("/bin/basename -s .out '" + fullfilepath + "'")
# ion_name = runbash("echo '" + filename + "' | cut -d'-' -f2")
# base_name = runbash("echo '" + filename + "' | cut -d'-' -f1")
# electrostatics = runbash("cat '" + i + "' | grep 'Electrostatics ' -m1 | tr -s ' ' | cut -d' ' -f7")
# exchange = runbash("cat '" + i + "' | grep 'Exchange ' -m1 | tr -s ' ' | cut -d' ' -f7")
# induction = runbash("cat '" + i + "' | grep 'Induction ' -m1 | tr -s ' ' | cut -d' ' -f7")
# dispersion = runbash("cat '" + i + "' | grep 'Dispersion ' -m1 | tr -s ' ' | cut -d' ' -f7")
# total = runbash("cat '" + i + "' | grep ' Total SAPT2+ ' -m1 | tr -s ' ' | cut -d' ' -f8")
# charge_transfer = runbash("cat '" + i + "' | grep 'SAPT Charge Transfer' | tail -n 1 | tr -s ' ' | cut -d' ' -f9")
# newdata = {"sort" : sort[base_name], "base_name" : base_name, "ion_name" : ion_name, "electrostatics" : electrostatics, "exchange" : exchange, "induction" : induction, "dispersion" : dispersion, "total" : total, "charge_transfer" : charge_transfer, "dhomo" : 0, "dlumo" : 0, "plotsize" : 0}
# datapairs = datapairs.append(newdata, ignore_index=True)
# # Drops any frames that have a "None" value in charge_trasfer (the last SAPT calculation)
# datapairs = datapairs[datapairs['charge_transfer'] != None]
# # Sorts the frames byt their name so that they are plotted alphabetically
# datapairs = datapairs.sort_values(["ion_name", "sort"])
# # Sets the dataframe column type to float
# datapairs.electrostatics = datapairs.electrostatics.astype(float)
# datapairs.exchange = datapairs.exchange.astype(float)
# datapairs.induction = datapairs.induction.astype(float)
# datapairs.dispersion = datapairs.dispersion.astype(float)
# datapairs.total = datapairs.total.astype(float)
# datapairs.charge_transfer = datapairs.charge_transfer.astype(float)
# datapairs.dhomo = datapairs.dhomo.astype(float)
# datapairs.dlumo = datapairs.dlumo.astype(float)
# datapairs.bhomo = datapairs.bhomo.astype(float)
# datapairs.blumo = datapairs.blumo.astype(float)
# datapairs.plotsize = datapairs.plotsize.astype(int)
# datapairs = pd.merge(datapairs, homodatapairs, on=['base_name', 'ion_name'])
# Eh2eV(datapairs, ["lumo", "homo", "dhomo", "dlumo", "bhomo", "blumo"])
# Eh2eV(lonepairs, ["lumo", "homo"])
# warnings.filterwarnings("ignore")
# for i in range(len(datapairs.base_name)):
# for j in range(len(lonepairs.homo)):
# if datapairs.base_name[i] == lonepairs.base_name[j]:
# datapairs.dhomo[i] = datapairs.homo[i] - lonepairs.homo[j]
# datapairs.dlumo[i] = datapairs.lumo[i] - lonepairs.lumo[j]
# datapairs.bhomo[i] = lonepairs.homo[j]
# datapairs.blumo[i] = lonepairs.lumo[j]
# for i in range(len(datapairs.base_name)):
# if math.isnan(datapairs.total[i]) == False:
# datapairs.plotsize[i] = abs(int(datapairs.total[i]))
# data.to_csv(r'./IonsData/pd.data.csv')
# datapairs.to_csv(r'./IonsData/pd.datapairs.csv')
data = pd.read_csv("./IonsData/pd.data.csv")
datapairs = pd.read_csv("./IonsData/pd.datapairs.csv")
data = data.sort_values(["ion_name", "sort"])
datapairs.electrostatics = datapairs.electrostatics.astype(float)
datapairs.exchange = datapairs.exchange.astype(float)
datapairs.induction = datapairs.induction.astype(float)
datapairs.dispersion = datapairs.dispersion.astype(float)
datapairs.total = datapairs.total.astype(float)
datapairs.charge_transfer = datapairs.charge_transfer.astype(float)
datapairs.dhomo = datapairs.dhomo.astype(float)
datapairs.dlumo = datapairs.dlumo.astype(float)
datapairs.bhomo = datapairs.bhomo.astype(float)
datapairs.blumo = datapairs.blumo.astype(float)
datapairs.plotsize = datapairs.plotsize.astype(int)
data.electrostatics = data.electrostatics.astype(float)
data.exchange = data.exchange.astype(float)
data.induction = data.induction.astype(float)
data.dispersion = data.dispersion.astype(float)
data.total = data.total.astype(float)
data.charge_transfer = data.charge_transfer.astype(float)
data.dhomo = data.dhomo.astype(float)
data.dlumo = data.dlumo.astype(float)
data.bhomo = data.bhomo.astype(float)
data.blumo = data.blumo.astype(float)
data.plotsize = data.plotsize.astype(int)
def plothomo(df, ionclass):
if ionclass == "metal":
widthval = 500
else:
widthval = 1000
fig = px.bar(df[df.ion_class == ionclass],
x="base_name",
y="homo",
facet_col='ion_name',
labels={'homo':'Energy (eV)', 'base_name':""},
height=400,
width=widthval,
barmode="overlay",
opacity=1,
)
for i in range(len(fig.data)):
fig.data[i].name="HOMO"
fig2 = px.bar(df[df.ion_class == ionclass],
x="base_name",
y="dhomo",
facet_col='ion_name',
labels={'homo':'Energy (eV)', 'base_name':""},
height=400,
width=widthval,
barmode="overlay",
opacity=1,
color_discrete_map={"dhomo": "rgb(255,119,119)"}
)
for i in range(len(fig2.data)):
fig2.data[i].name="ΔHOMO"
fig.add_trace(fig2.data[0])
fig.add_trace(fig2.data[1])
if ionclass != "metal":
fig.add_trace(fig2.data[2])
fig.add_trace(fig2.data[3])
fig.add_trace(fig2.data[4])
fig.add_trace(fig2.data[5])
if ionclass != "cation":
fig.add_trace(fig2.data[6])
fig.data[1].showlegend = True
fig.data[-1].showlegend = True
fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
return(fig.show())
for ionclass in classnames:
plothomo(data, ionclass)
Excluding Benzene-Ni
def saptdf(df):
newdf = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'energy', 'energytype'])
for i in list(range(len(df))):
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.total[i],"energytype" : "total"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.electrostatics[i],"energytype" : "electrostatics"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.exchange[i],"energytype" : "exchange"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.induction[i],"energytype" : "induction"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.dispersion[i],"energytype" : "dispersion"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.charge_transfer[i],"energytype" : "charge_transfer"},ignore_index=True)
newdf.sort_values(["ion_name", "sort"])
return(newdf)
def plotsapt(df, ionclass):
dfint = saptdf(df)
if ionclass == "metal":
widthval = 550
else:
widthval = 1000
fig = px.scatter(dfint[dfint.ion_class == ionclass],
x="base_name",
y="energy",
facet_col='ion_name',
labels={'energy':'Energy (KJ/mol)', 'base_name':""},
height=400,
width=widthval,
# barmode="overlay",
# opacity=0.5,
color="energytype",
)
fig.layout.showlegend = True
fig.update_traces(opacity=0.7, marker=dict(size=10, line=dict(color='black', width=1)))
fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
fig.for_each_trace(lambda t: t.update(name=t.name.replace("energytype=","")))
return(fig.show())
for ionclass in classnames:
plotsapt(data, ionclass)
trendlist = ["ion_name", "ion_class", "base_name"]
Excluding Pd$^{2+}$, Ni$^{2+}$ and hard ions(Cl$^-$, Na$^+$ and Li$^+$)
dplusel is dispersion + electrostatics
# Feature Plot
datacorr = data
datacorr['dplusel'] = np.add(datacorr['electrostatics'], datacorr['dispersion'])
# Correlation plot
fig = px.scatter_matrix(datacorr[datacorr["ion_class"] != "metal"][datacorr["ion_name"] != "Na"][datacorr["ion_name"] != "Li"][datacorr["ion_name"] != "Cl"],
dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
color="ion_name",
hover_data=['ion_name', 'base_name'],
width=1000,
height=1000)
fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()
Excluding Pd$^{2+}$, Ni$^{2+}$ and hard ions(Cl$^-$, Na$^+$ and Li$^+$)
dplusel is dispersion + electrostatics
datacorr = data
datacorr['dplusel'] = np.add(datacorr['electrostatics'], datacorr['dispersion'])
datacorr['dpluselplusex'] = np.add(datacorr['dplusel'], datacorr['exchange'])
# Correlation plot
fig = px.scatter_matrix(datacorr[datacorr["ion_class"] != "metal"][datacorr["ion_name"] != "Na"][datacorr["ion_name"] != "Li"][datacorr["ion_name"] != "Cl"],
dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
color="ion_class",
hover_data=['ion_name', 'base_name'],
width=1000,
height=1000)
fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()
def plothomopairs(df):
fig = px.bar(df,
x="base_name",
y="homo",
facet_col='ion_name',
labels={'homo':'Energy (eV)', 'base_name':""},
width=600,
height=400,
barmode="overlay",
opacity=1,
)
for i in range(len(fig.data)):
fig.data[i].name="HOMO"
fig2 = px.bar(df,
x="base_name",
y="dhomo",
facet_col='ion_name',
labels={'homo':'Energy (eV)', 'base_name':""},
height=400,
barmode="overlay",
opacity=1,
color_discrete_map={"dhomo": "rgb(31,120,180)"}
)
for i in range(len(fig2.data)):
fig2.data[i].name="ΔHOMO"
fig.add_trace(fig2.data[0])
fig.add_trace(fig2.data[1])
fig.data[1].showlegend = True
fig.data[-1].showlegend = True
fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
return(fig.show())
plothomopairs(datapairs)
def saptdfpairs(df):
newdf = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'energy', 'energytype'])
for i in list(range(len(df))):
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.total[i],"energytype" : "total"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.electrostatics[i],"energytype" : "electrostatics"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.exchange[i],"energytype" : "exchange"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.induction[i],"energytype" : "induction"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.dispersion[i],"energytype" : "dispersion"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.charge_transfer[i],"energytype" : "charge_transfer"},ignore_index=True)
newdf.sort_values(["ion_name", "sort"])
return(newdf)
def plotsaptpairs(df):
dfint = saptdfpairs(df)
if ionclass == "metal":
widthval = 550
else:
widthval = 600
fig = px.scatter(dfint,
x="base_name",
y="energy",
facet_col='ion_name',
labels={'homo':'Energy (eV)', 'base_name':""},
width=600,
height=400,
# width=widthval,
# barmode="overlay",
color="energytype",
)
fig.layout.showlegend = True
fig.update_traces(opacity=0.7, marker=dict(size=12, line=dict(color='black', width=1)))
fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
fig.for_each_trace(lambda t: t.update(name=t.name.replace("energytype=","")))
return(fig.show())
plotsaptpairs(datapairs)
databoth = pd.DataFrame(columns=["base_name", "ion_name", "electrostatics", "exchange", "induction", "dispersion", "total", "charge_transfer", "dhomo", "sorp"])
for i in range(len(list(data.base_name))):
newdata = {"base_name":data.base_name[i], "ion_name":data.ion_name[i], "electrostatics":data.electrostatics[i], "exchange":data.exchange[i], "induction":data.induction[i], "dispersion":data.dispersion[i], "total":data.total[i], "charge_transfer":data.charge_transfer[i], "dhomo":data.dhomo[i], "sorp":"s"}
databoth = databoth.append(newdata, ignore_index=True)
for i in range(len(list(datapairs.base_name))):
newdata = {"base_name":datapairs.base_name[i], "ion_name":datapairs.ion_name[i], "electrostatics":datapairs.electrostatics[i], "exchange":datapairs.exchange[i], "induction":datapairs.induction[i], "dispersion":datapairs.dispersion[i], "total":datapairs.total[i], "charge_transfer":datapairs.charge_transfer[i], "dhomo":datapairs.dhomo[i], "sorp":"p"}
databoth = databoth.append(newdata, ignore_index=True)
databoth = databoth[databoth["ion_name"] != "Ni"][databoth["ion_name"] != "Pd"][databoth["ion_name"] != "Na"][databoth["ion_name"] != "Cl"][databoth["ion_name"] != "Li"][databoth["ion_name"] != "Na"]
databoth['dplusel'] = np.add(databoth['electrostatics'], databoth['dispersion'])
databoth['dpluselplusex'] = np.add(databoth['dplusel'], databoth['exchange'])
fig = px.scatter_matrix(databoth,
dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
color="sorp",
hover_data=['ion_name', 'base_name'],
width=1000,
height=1000)
fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()
Except Li (Na is missing for some reason) the electrostatics follows the total interaction energy. Could you please extract the gradient from the electrostatics-total graph. I believe that Induction + dispersion + exchange might be summing to within a few kJ mol-1.... Ultimately, electrostatics is the force that drives this interaction... Li definitely makes the bond covalent and hence such a huge difference between electrostatics and total energy.
Total - Electrostatics¶
Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$¶
data['dplusel'] = np.add(data['electrostatics'], data['dispersion'])
data['dpluselplusex'] = np.add(data['dplusel'], data['exchange'])
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
y="total",
x="electrostatics",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("total = " + str(results.px_fit_results.iloc[0].params[1]) +"* Electrostatics " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.9030553646414368 total = 0.9856406966724351* Electrostatics -6.433560778791372
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
y="total",
x="electrostatics",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Total = " + str(results.px_fit_results.iloc[0].params[1]) +" * Electrostatics + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.8904670996931319 Total = 0.8736926799703468 * Electrostatics + -7.35623462676404
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
y="total",
x="dplusel",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("total = " + str(results.px_fit_results.iloc[0].params[1]) +"* (Electrostatics + Dispersion) " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.6289566526512076 total = 0.743313858961744* (Electrostatics + Dispersion) 3.291513720534722
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
y="total",
x="dplusel",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("total = " + str(results.px_fit_results.iloc[0].params[1]) +"* (Electrostatics + Dispersion) " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.7681970967915128 total = 0.6042561976435759* (Electrostatics + Dispersion) 1.9471283622225952
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
y="total",
x="dpluselplusex",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Total = " + str(results.px_fit_results.iloc[0].params[1]) +" * (Electrostatics + Dispersion + Exchange) + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.7918257575866938 Total = 1.7648848320532988 * (Electrostatics + Dispersion + Exchange) + -5.364755450287623
fig = px.scatter(data[datacorr.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
y="total",
x="dpluselplusex",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Total = " + str(results.px_fit_results.iloc[0].params[1]) +" * (Electrostatics + Dispersion + Exchange) + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.9424319306598973 Total = 1.4124998291175728 * (Electrostatics + Dispersion + Exchange) + -6.249230648555187
Electrostatics and exchange also correlate with each other, and so do a) dispersion and exchange and b) induction and electrostatics. Extracting their gradients would be good too. Surprisingly, dispersion correlates with exchange without the need to add electrostatics! Although the latter does improve the correlation.... The dispersion-exchange correlation could be compared to Anh's data.... I think this is the reason why her DFT works so well...
Electrostatics - Exchange¶
Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$¶
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
y="electrostatics",
x="exchange",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Electrostatics = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.5118609630731161 Electrostatics = -1.063969275341984 * Exchange + 9.512563076218635
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
y="electrostatics",
x="exchange",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Electrostatics = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.5316140446902322 Electrostatics = -0.8780190254123053 * Exchange + 8.361633853052801
fig = px.scatter(datacorr[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
y="dispersion",
x="exchange",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Dispersion = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.37000294355481467 Dispersion = -0.4877974627087508 * Exchange + -6.003526161695268
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
y="dispersion",
x="exchange",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Dispersion = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.7464544651613487 Dispersion = -0.6360920365032324 * Exchange + -4.788495179125986
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
y="induction",
x="electrostatics",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Induction = " + str(results.px_fit_results.iloc[0].params[1]) +" * Electrostatics + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.6506407144336758 Induction = 0.4998111176981335 * Electrostatics + -4.337012515738229
fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
y="induction",
x="electrostatics",
# color="sorp",
height=400,
width=600,
hover_data=['ion_name', 'base_name'],
trendline="ols",
# size='plotsize'
)
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Induction = " + str(results.px_fit_results.iloc[0].params[1]) +" * Electrostatics + " + str(results.px_fit_results.iloc[0].params[0]))
R^2 = 0.5906389608317781 Induction = 0.26427779326417417 * Electrostatics + -6.42340309311744
I would like to know where the ion pairs are located on 3.1.1 plots. Could you please include only those of the lowest energy? Would you be able to make SAPT energy graphs combining the ion pairs and the results for individual ions included in these ion pairs?
databoth2 = pd.DataFrame(columns=["base_name", "ion_name", "electrostatics", "exchange", "induction", "dispersion", "total", "charge_transfer", "dhomo", "sorp"])
for i in range(len(list(data.base_name))):
newdata = {"base_name":data.base_name[i], "ion_name":data.ion_name[i], "electrostatics":data.electrostatics[i], "exchange":data.exchange[i], "induction":data.induction[i], "dispersion":data.dispersion[i], "total":data.total[i], "charge_transfer":data.charge_transfer[i], "dhomo":data.dhomo[i], "sorp":"s"}
databoth2 = databoth2.append(newdata, ignore_index=True)
for i in range(len(list(datapairs.base_name))):
newdata = {"base_name":datapairs.base_name[i], "ion_name":datapairs.ion_name[i], "electrostatics":datapairs.electrostatics[i], "exchange":datapairs.exchange[i], "induction":datapairs.induction[i], "dispersion":datapairs.dispersion[i], "total":datapairs.total[i], "charge_transfer":datapairs.charge_transfer[i], "dhomo":datapairs.dhomo[i], "sorp":"p"}
databoth2 = databoth2.append(newdata, ignore_index=True)
databoth2 = databoth2[databoth2.ion_name != "BF2(CF3)2"][databoth2.ion_name != "C(CN)3"][databoth2.ion_name != "N(SOOCF3)2"][databoth2.ion_name != "N(NC)2"][databoth2.ion_name != "C1MPyr"][databoth2.ion_name != "NMe4"][databoth2.ion_name != "PMe4"][databoth2.ion_name != "Na"][databoth2.ion_name != "Pd"][databoth2.ion_name != "Ni"][databoth2.ion_name != "tempo"]
databoth2['dplusel'] = np.add(databoth2['electrostatics'], databoth2['dispersion'])
databoth2['dpluselplusex'] = np.add(databoth2['dplusel'], databoth2['exchange'])
databoth2
fig = px.scatter_matrix(databoth2,
dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
color="sorp",
hover_data=['ion_name', 'base_name'],
width=1000,
height=1000)
fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()
Ni and Pd form covalent bonds by the look of it. In both cases it is the induction that dominates the interaction (hence, the symmetry is crucial). Electrostatics is stronger in Ni but is compensated by high exchange... So induction is the key whereas in ILs pure electrostatics is the key! It must come down to the orbital symmetry I guess. This is very cool though as clearly transition metals favour certain sites whereas the problem of ions is that they lead to non-directional interaction... So, the next goal would be to identify ions that have larger induction over electrostatics (if it is even possible).
This problem seems to be dependent on the molecule with which the ion is interacting and unfortunately we can only really say that BF4$^+$ and Li$^-$ have the highest number of incidence where induction is greater than electrostatics.
def saptdf(df):
newdf = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'energy', 'energytype'])
for i in list(range(len(df))):
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.electrostatics[i],"energytype" : "electrostatics"},ignore_index=True)
newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.induction[i],"energytype" : "induction"},ignore_index=True)
newdf.sort_values(["ion_name", "sort"])
return(newdf)
def plotsapt(df, ionclass):
dfint = saptdf(df)
if ionclass == "metal":
widthval = 550
else:
widthval = 1000
fig = px.scatter(dfint[dfint.ion_class == ionclass],
x="base_name",
y="energy",
facet_col='ion_name',
labels={'energy':'Energy (KJ/mol)', 'base_name':""},
height=400,
width=widthval,
# barmode="overlay",
# opacity=0.5,
color="energytype",
)
fig.layout.showlegend = True
fig.update_traces(opacity=0.7, marker=dict(size=10, line=dict(color='black', width=1)))
fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
fig.for_each_trace(lambda t: t.update(name=t.name.replace("energytype=","")))
return(fig.show())
for ionclass in classnames:
plotsapt(data, ionclass)